#include <dftHardSphere.h>

// Note: '2' as in 'rho2' or 'rho2Hat' or 'WHat2' refers to vectors living on reflected DAs 
// (twice the size in each dimension of the "physical" DAs.
// Likewise for the actual DAs.


static PetscErrorCode CalculateWHat2_1d(DFTHardSphere *dft);
static PetscErrorCode CreateFFT2(DFTHardSphere* dft);
static PetscErrorCode DestroyFFT2(DFTHardSphere* dft);



#undef __FUNCT__
#define __FUNCT__ "DFTHardSphereCreate"
PetscErrorCode DFTHardSphereCreate(DFTHardSphere *dft) {

  PetscErrorCode ierr;

  PetscFunctionBegin;
  
  PetscTruth flag;
  dft->check_pf = PETSC_TRUE;
  ierr = PetscOptionsGetTruth("dft_hs_", "-check_pf", &(dft->check_pf), &flag); CHKERRQ(ierr);

  // Only 1d for now
  if(dft->dim == 1) {
    dft->dims2[0] = 2.0*dft->dims[0];
    ierr = DACreate1d(dft->comm, DA_NONPERIODIC, dft->dims2[0], NUM_SPECIES, 0, PETSC_NULL, &(dft->rhoDA2)); CHKERRQ(ierr);
    ierr = DACreate1d(dft->comm, DA_NONPERIODIC, dft->dims[0], 1, 0, PETSC_NULL, &(dft->singleDA)); CHKERRQ(ierr);
    ierr = DACreate1d(dft->comm, DA_NONPERIODIC, dft->dims2[0], 1, 0, PETSC_NULL, &(dft->singleDA2)); CHKERRQ(ierr);
    dft->numAlpha = 6;
    ierr = DACreate1d(dft->comm, DA_NONPERIODIC, dft->dims[0], dft->numAlpha, 0, PETSC_NULL, &(dft->weightsDA)); CHKERRQ(ierr);
    ierr = DACreate1d(dft->comm, DA_NONPERIODIC, dft->dims2[0], dft->numAlpha, 0, PETSC_NULL, &(dft->weightsDA2)); CHKERRQ(ierr);
  }
  else {
    SETERRQ1(PETSC_ERR_SUP, "Unsupported domain dimension: %d", dft->dim);
  }
  ierr = CreateFFT2(dft); CHKERRQ(ierr);

  /* Create auxiliary Vecs */
  /** Rho DAs **/
  ierr = DAGetGlobalVector(dft->rhoDA2, &(dft->rho2)); CHKERRQ(ierr);
  ierr = DAGetGlobalVector(dft->rhoDA2, &(dft->mu2));  CHKERRQ(ierr);

  /** Weights DAs **/
  ierr = DAGetGlobalVector(dft->weightsDA2, &(dft->n2)); CHKERRQ(ierr);
  ierr = DAGetGlobalVector(dft->weightsDA2, &(dft->PhiDer2)); CHKERRQ(ierr);
  ierr = DAGetGlobalVector(dft->weightsDA2, &(dft->PhiDer2Hat)); CHKERRQ(ierr);

  /** Single DAs **/
  // z is a purely auxiliary accumulator Vec
  ierr = DAGetGlobalVector(dft->singleDA2, &(dft->z2)); CHKERRQ(ierr);

  /* Create auxiliary split Vec sets */
  ierr = PetscMalloc5(dft->numSpecies, Vec, &(dft->rho2BySpecies), 
                      dft->numSpecies, Vec, &(dft->rho2HatBySpecies), 
                      dft->numSpecies, Vec, &(dft->mu2BySpecies),
                      dft->numSpecies, Vec, &(dft->mu2HatBySpecies),
                      dft->numSpecies, Vec, &(dft->WHat2BySpecies));  
  CHKERRQ(ierr);
  for(PetscInt sp = 0; sp < dft->numSpecies; ++sp) { 
    ierr = DAGetGlobalVector(dft->singleDA2,  &(dft->rho2BySpecies[sp]));    CHKERRQ(ierr);
    ierr = DAGetGlobalVector(dft->singleDA2,  &(dft->rho2HatBySpecies[sp])); CHKERRQ(ierr);
    ierr = DAGetGlobalVector(dft->singleDA2,  &(dft->mu2BySpecies[sp]));     CHKERRQ(ierr);     
    ierr = DAGetGlobalVector(dft->singleDA2,  &(dft->mu2HatBySpecies[sp]));  CHKERRQ(ierr);     
    ierr = DAGetGlobalVector(dft->weightsDA2, &dft->WHat2BySpecies[sp]);     CHKERRQ(ierr);
  } 
  ierr = PetscMalloc5(dft->numAlpha, Vec, &(dft->WHat2SpeciesByAlpha),
                      dft->numAlpha, Vec, &(dft->n2ByAlpha),
                      dft->numAlpha, Vec, &(dft->n2HatByAlpha),
                      dft->numAlpha, Vec, &(dft->PhiDer2ByAlpha), 
                      dft->numAlpha, Vec, &(dft->PhiDer2HatByAlpha));
  CHKERRQ(ierr);
  for(PetscInt alpha = 0; alpha < dft->numAlpha; ++alpha) { 
    ierr = DAGetGlobalVector(dft->singleDA2, &(dft->WHat2SpeciesByAlpha[alpha])); CHKERRQ(ierr);
    ierr = DAGetGlobalVector(dft->singleDA2, &(dft->n2ByAlpha[alpha]));           CHKERRQ(ierr);     
    ierr = DAGetGlobalVector(dft->singleDA2, &(dft->n2HatByAlpha[alpha]));        CHKERRQ(ierr);     
    ierr = DAGetGlobalVector(dft->singleDA2, &(dft->PhiDer2ByAlpha[alpha]));      CHKERRQ(ierr);
    ierr = DAGetGlobalVector(dft->singleDA2, &(dft->PhiDer2HatByAlpha[alpha]));   CHKERRQ(ierr);
  } 
  ierr = CalculateWHat2_1d(dft); CHKERRQ(ierr);
  PetscFunctionReturn(0);

}//DFTHardSphereCreate()

#undef __FUNCT__
#define __FUNCT__ "DFTHardSphereDestroy"
PetscErrorCode DFTHardSphereDestroy(DFTHardSphere *dft)
{

  PetscErrorCode ierr;
  PetscFunctionBegin;

  /* Destroy auxiliary vectors */
  /** Rho DAs **/
  ierr = DARestoreGlobalVector(dft->rhoDA2, &(dft->rho2)); CHKERRQ(ierr);
  ierr = DARestoreGlobalVector(dft->rhoDA2, &(dft->mu2));  CHKERRQ(ierr);

  /** Weights DAs **/
  ierr = DARestoreGlobalVector(dft->weightsDA2, &(dft->PhiDer2Hat)); CHKERRQ(ierr);
  ierr = DARestoreGlobalVector(dft->weightsDA2,   &(dft->PhiDer2)); CHKERRQ(ierr);
  ierr = DARestoreGlobalVector(dft->weightsDA2,   &(dft->n2));      CHKERRQ(ierr);


  /** Single DAs **/
  ierr = DARestoreGlobalVector(dft->singleDA2, &(dft->z2)); CHKERRQ(ierr);

  /* Destroy auxiliary split Vec sets */
  for(PetscInt sp = 0; sp < dft->numSpecies; ++sp) { 
    ierr = DARestoreGlobalVector(dft->singleDA2,  &(dft->rho2BySpecies[sp]));     CHKERRQ(ierr);
    ierr = DARestoreGlobalVector(dft->singleDA2,  &(dft->rho2HatBySpecies[sp]));  CHKERRQ(ierr);
    ierr = DARestoreGlobalVector(dft->singleDA2,  &(dft->mu2BySpecies[sp]));      CHKERRQ(ierr);     
    ierr = DARestoreGlobalVector(dft->singleDA2,  &(dft->mu2HatBySpecies[sp]));   CHKERRQ(ierr);     
    ierr = DARestoreGlobalVector(dft->weightsDA2, &(dft->WHat2BySpecies[sp]));    CHKERRQ(ierr);     
  } 
  ierr = PetscFree5((dft->rho2BySpecies), 
                    (dft->rho2HatBySpecies), 
                    (dft->mu2BySpecies),
                    (dft->mu2HatBySpecies),
                    (dft->WHat2BySpecies));  
  CHKERRQ(ierr);
  for(PetscInt alpha = 0; alpha < dft->numAlpha; ++alpha) { 
    ierr = DARestoreGlobalVector(dft->singleDA2, &(dft->WHat2SpeciesByAlpha[alpha])); CHKERRQ(ierr);
    ierr = DARestoreGlobalVector(dft->singleDA2, &(dft->n2ByAlpha[alpha]));           CHKERRQ(ierr);
    ierr = DARestoreGlobalVector(dft->singleDA2, &(dft->n2HatByAlpha[alpha]));        CHKERRQ(ierr);
    ierr = DARestoreGlobalVector(dft->singleDA2, &(dft->PhiDer2ByAlpha[alpha]));      CHKERRQ(ierr);
    ierr = DARestoreGlobalVector(dft->singleDA2, &(dft->PhiDer2HatByAlpha[alpha]));   CHKERRQ(ierr);
  } 
  ierr = PetscFree5((dft->WHat2SpeciesByAlpha),
                    (dft->n2ByAlpha),
                    (dft->n2HatByAlpha),
                    (dft->PhiDer2ByAlpha), 
                    (dft->PhiDer2HatByAlpha));

  CHKERRQ(ierr);


  ierr = DestroyFFT2(dft); CHKERRQ(ierr);

  /* Destroy DAs */
  ierr = DADestroy(dft->weightsDA); CHKERRQ(ierr);
  ierr = DADestroy(dft->weightsDA2); CHKERRQ(ierr);
  ierr = DADestroy(dft->singleDA); CHKERRQ(ierr);
  ierr = DADestroy(dft->singleDA2); CHKERRQ(ierr);
  ierr = DADestroy(dft->rhoDA2); CHKERRQ(ierr);
  PetscFunctionReturn(0);

}//DFTHardSphereDestroy()

// sc is a name I made up for the function below (essentially, sinc) to remove its apparent singularity
inline PetscScalar sc(const PetscReal k, const PetscReal R) {
  PetscReal u = k*R;
  return (u== 0.0)?1.0:sin(u)/u;
}// sc()

// cs is a name I made up for the function below to remove its apparent singularity
inline PetscScalar cs(const PetscReal k, const PetscReal R) {
  // use expansions for sin and cos at 0 to figure out the leading order term
  PetscReal kR = k*R;
  return (kR!=0.0)?(2.0*R/(k*k))*(cos(kR)/(k*k)-sc(k,R)):-((2.0*R*R*R)/3.0);
}// cs()

PetscScalar W_0(const DFTHardSphere* dft, const PetscReal xx[], const int species) {
  return 1.0/(2.0*dft->R[species]);
}

PetscScalar WHat_0(const DFTHardSphere* dft, const PetscReal kk[], const int species) {
  return (2.0*sc(kk[0],dft->R[species]))/(2.0*sqrt(2*M_PI));
}

PetscScalar W_1(const DFTHardSphere* dft, const PetscReal xx[], const int species) {
  return 0.5;
}

PetscScalar WHat_1(const DFTHardSphere* dft, const PetscReal kk[], const int species) {
  return 0.5*(2.0*dft->R[species]*sc(kk[0],dft->R[species]))/(sqrt(2*M_PI));
}

PetscScalar W_2(const DFTHardSphere* dft, const PetscReal xx[], const int species) {
  return 2.0*M_PI*dft->R[species];
}

PetscScalar WHat_2(const DFTHardSphere* dft, const PetscReal kk[], const int species) {
  return 2.0*M_PI*dft->R[species]*(2.0*dft->R[species]*sc(kk[0],dft->R[species]))/(sqrt(2*M_PI));
}

PetscScalar W_3(const DFTHardSphere* dft, const PetscReal xx[], const int species) {
  return M_PI*(dft->R[species]*dft->R[species] - xx[0]*xx[0]);
}

PetscScalar WHat_3(const DFTHardSphere* dft, const PetscReal kk[], const int species) {
  PetscReal R = dft->R[species], k = kk[0];
  return M_PI*((R*R*R*2.0*sc(k,R))/(sqrt(2*M_PI)) - (2.0/sqrt(2.0*M_PI))*((R*R*R*sc(k,R))+cs(k,R)));
}

PetscScalar W_V1(const DFTHardSphere* dft, const PetscReal xx[], const int species) {
  return 0.5*xx[0]/dft->R[species];
}

PetscScalar WHat_V1(const DFTHardSphere* dft, const PetscReal kk[], const int species) {
  PetscScalar i(0.0,1.0);
  PetscReal k = kk[0], R = dft->R[species];
  return (0.5/R)*(i/sqrt(2.0*M_PI))*(k*cs(k,R));
}

PetscScalar W_V2(const DFTHardSphere* dft, const PetscReal xx[], const int species) {
  return 2.0*M_PI*xx[0];
}

PetscScalar WHat_V2(const DFTHardSphere* dft, const PetscReal kk[], const int species) {
  PetscScalar i(0.0,1.0);
  PetscReal k = kk[0], R = dft->R[species];
  return 2.0*M_PI*(i/sqrt(2.0*M_PI))*(k*cs(k,R));
}

#undef __FUNCT__
#define __FUNCT__ "CalculateWHat2_1d"
PetscErrorCode CalculateWHat2_1d(DFTHardSphere *dft)
{
  PetscScalar     (*WHatFunc[6])(const DFTHardSphere *, const PetscReal [], const int);
  PetscErrorCode    ierr;

  PetscFunctionBegin;
  if (dft->numAlpha != 6) {
    SETERRQ2(PETSC_ERR_SUP, "Invalid bound for alpha in 1d: %d should be %d", dft->numAlpha, 6);
  }
  WHatFunc[0] = WHat_0;
  WHatFunc[1] = WHat_1;
  WHatFunc[2] = WHat_2;
  WHatFunc[3] = WHat_3;
  WHatFunc[4] = WHat_V1;
  WHatFunc[5] = WHat_V2;
  
  for(PetscInt sp = 0; sp < dft->numSpecies; ++sp) {
    Weights *WHat2BySpeciesArray;
    ierr = DAVecGetArray(dft->weightsDA2, dft->WHat2BySpecies[sp], &WHat2BySpeciesArray); CHKERRQ(ierr);
    /* 
       I claim this is a better way of localizing this computation: disregard the actual physical coordinate,
       and go solely by the local wavenumber.  That way we might be able to use local Fourier exapansions
       in the future
    */
    /* For now I'll use global coordinates and domain sizes.  Or maybe dft->dims should be global? */
    PetscReal kk[1];
    for(PetscInt i = 0; i < dft->dims2[0]; ++i) {
      if(i <= dft->dims[0]) { // should be if(i <= dft->dims2[0]), but dims2[0] = 2.0*dims[0] by design 
        kk[0] = 2.0*M_PI*i/(2.0*dft->L[0]);  // can reduce this fraction by 2.0, but leave as is for clarity
      } 
      else {
        kk[0] = -2.0*M_PI*(dft->dims2[0]-i)/(2.0*dft->L[0]); // can reduce this fraction by 2.0, but leave as is for clarity
      }
      for(PetscInt alpha = 0; alpha < dft->numAlpha; ++alpha) {
        WHat2BySpeciesArray[i].w[alpha] = WHatFunc[alpha](dft, kk, sp);
        if(dft->debug > DFT_DEBUG_HARD_SPHERE_WHAT2_DETAIL) {
          ierr = PetscPrintf(PETSC_COMM_WORLD, "%s: sp = %d, i = %d, alpha = %d, kk[0] = %g, WHat2BySpeciesArray[%d].w[%d] = %g+i*%g: \n", 
                             __FUNCT__, sp, i, alpha, kk[0], i, alpha, 
                             PetscRealPart(WHat2BySpeciesArray[i].w[alpha]), 
                             PetscImaginaryPart(WHat2BySpeciesArray[i].w[alpha])); 
          CHKERRQ(ierr);
        }
      }//alpha-loop
    }// ii-loop
    ierr = DAVecRestoreArray(dft->weightsDA2, dft->WHat2BySpecies[sp], &WHat2BySpeciesArray); CHKERRQ(ierr);
    if(dft->debug > DFT_DEBUG_HARD_SPHERE_WHAT2) {
      ierr = PetscPrintf(PETSC_COMM_WORLD, "%s: WHat2BySpecies[%d]: \n", __FUNCT__, sp); CHKERRQ(ierr);
      ierr = VecView(dft->WHat2BySpecies[sp], PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(ierr);
    }
  }// sp-loop
  PetscFunctionReturn(0);
}// CalculateWHat2_1d()




#undef __FUNCT__
#define __FUNCT__ "CreateFFT2"
PetscErrorCode CreateFFT2(DFTHardSphere *dft)
{
  PetscErrorCode ierr;
  PetscFunctionBegin;
  // Note FFTW wants dimensions in reverse order (since we have Fortran storage)
  PetscInt dims[3];
  if(dft->dim == 1) {
    dims[0] = 2*dft->dims[0];
    ierr = MatCreateSeqFFTW(PETSC_COMM_WORLD, dft->dim, dims, &dft->F2);CHKERRQ(ierr);
  }
  else if(dft->dim == 3) {
    PetscInt dims[3];
    dims[0] = 2*dft->dims[2];
    dims[1] = 2*dft->dims[1];
    dims[2] = 2*dft->dims[0];
    ierr = MatCreateSeqFFTW(PETSC_COMM_WORLD, dft->dim, dims, &dft->F2);CHKERRQ(ierr);
  }
  else {
    SETERRQ1(PETSC_ERR_SUP, "Unsupported domain dimension: %d", dft->dim);
  }  
  PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "DestroyFFT2"
PetscErrorCode DestroyFFT2(DFTHardSphere* dft)
{
  PetscErrorCode ierr;
  PetscFunctionBegin;
  ierr = MatDestroy(dft->F2);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}


#undef __FUNCT__
#define __FUNCT__ "Calculate_n_1d"
PetscErrorCode Calculate_n_1d(DFTHardSphere *dft, Vec rho2, Vec n2) {
  /* Note: rho2 lives over rhoDA2, n2 lives over weightsDA2 */

  PetscErrorCode ierr;
  PetscFunctionBegin; 

  ierr = VecStrideGatherAll(rho2, dft->rho2BySpecies, INSERT_VALUES);CHKERRQ(ierr); 
  for(PetscInt sp = 0; sp < dft->numSpecies; ++sp) {
    ierr = MatMult(dft->F2, dft->rho2BySpecies[sp], dft->rho2HatBySpecies[sp]);CHKERRQ(ierr);
  }
  for(PetscInt alpha = 0; alpha < dft->numAlpha; ++alpha) {
    ierr = VecSet(dft->n2HatByAlpha[alpha],0.0); CHKERRQ(ierr);
  }
  for(PetscInt sp = 0; sp < dft->numSpecies; ++sp) {
    ierr = VecStrideGatherAll(dft->WHat2BySpecies[sp], dft->WHat2SpeciesByAlpha, INSERT_VALUES); CHKERRQ(ierr);
    for(PetscInt alpha = 0; alpha < dft->numAlpha; ++alpha) {
      // Note: z2 is a purely auxiliary accumulator vector.
      VecPointwiseMult(dft->z2, dft->WHat2SpeciesByAlpha[alpha], dft->rho2HatBySpecies[sp]); CHKERRQ(ierr);
      ierr = VecAXPY(dft->n2HatByAlpha[alpha], 1.0, dft->z2);CHKERRQ(ierr); 
      PetscReal scale = 1.0/dft->dims2[0];
      // Correct for sign difference in convolution of odd functions
      if ((alpha == 4) || (alpha == 5))  {
        ierr = VecScale(dft->n2HatByAlpha[alpha], -scale);CHKERRQ(ierr);
      } else {
        ierr = VecScale(dft->n2HatByAlpha[alpha], scale);CHKERRQ(ierr);
      }
    }// alpha-loop
  }// sp-loop
  for(PetscInt alpha = 0; alpha < dft->numAlpha; ++alpha) {
    /* n2 back to physical space, one alpha at a time */
    ierr = MatMultTranspose(dft->F2, dft->n2HatByAlpha[alpha], dft->n2ByAlpha[alpha]); CHKERRQ(ierr);
    if(dft->check_pf && alpha == 3){
      // Check for valid packing fraction
      DALocalInfo    info;
      PetscScalar *n3;
      ierr = DAGetLocalInfo(dft->singleDA2, &info);CHKERRQ(ierr);
      ierr = DAVecGetArray(dft->singleDA2, dft->n2ByAlpha[alpha], &n3);CHKERRQ(ierr);
      for(PetscInt i = info.xs; i < info.xs+info.xm; ++i) {
        if (PetscRealPart(n3[i]) >= 1.0) {
          SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE, "Invalid packing fraction: n3[%d]: %g", i, PetscRealPart(n3[i]));
        }
      ierr = DAVecRestoreArray(dft->singleDA2, dft->n2ByAlpha[alpha], &n3);CHKERRQ(ierr);
      }
    }
  }// alpha-loop
  ierr = VecStrideScatterAll(dft->n2ByAlpha, n2, INSERT_VALUES);CHKERRQ(ierr); 
  PetscFunctionReturn(0); 
}/* Calcuate_n_1d() */


PetscScalar PhiDer_0(const PetscScalar n[]) {
  return -log(1.0 - PetscRealPart(n[3]));
}

PetscScalar PhiDer_1(const PetscScalar n[]) {
  return PetscRealPart(n[2])/(1.0 - PetscRealPart(n[3]));
}

PetscScalar PhiDer_2(const PetscScalar n[]) {
  const PetscScalar ratio = ((PetscRealPart(n[2]) == 0.0) && (PetscRealPart(n[5]) == 0.0)) ? 0.0 : PetscSqr(n[5])/PetscSqr(n[2]);
  const PetscScalar term  = 1.0 - ratio;
  return PetscRealPart(n[1])/(1.0 - PetscRealPart(n[3])) + (PetscSqr(PetscRealPart(n[2]))/(8.0*M_PI*PetscSqr(1.0 - PetscRealPart(n[3]))))*term*term*term + 
                     (PetscSqr(PetscRealPart(n[5]))/(4.0*M_PI*PetscSqr(1.0 - PetscRealPart(n[3]))))*PetscSqr(term);
}

PetscScalar PhiDer_3(const PetscScalar n[]) {
  const PetscScalar ratio = ((PetscRealPart(n[2]) == 0.0) && (PetscRealPart(n[5]) == 0.0)) ? 0.0 : PetscSqr(PetscRealPart(n[5]))/PetscSqr(PetscRealPart(n[2]));
  const PetscScalar term  = 1.0 - ratio;
  return PetscRealPart(n[0])/(1.0 - PetscRealPart(n[3])) + 
    (PetscRealPart(n[1])*PetscRealPart(n[2]) - PetscRealPart(n[4])*PetscRealPart(n[5]))/PetscSqr(1.0 - PetscRealPart(n[3])) + 
    (PetscRealPart(n[2])*PetscRealPart(n[2])*PetscRealPart(n[2])/(12.0*M_PI*(1.0 - PetscRealPart(n[3]))*(1.0 - PetscRealPart(n[3]))*(1.0 - PetscRealPart(n[3]))))*term*term*term;
}

PetscScalar PhiDer_V1(const PetscScalar n[]) {
  return -PetscRealPart(n[5])/(1.0 - PetscRealPart(n[3]));
}

PetscScalar PhiDer_V2(const PetscScalar n[]) {
  const PetscScalar ratio = ((PetscRealPart(n[2]) == 0.0) && (PetscRealPart(n[5]) == 0.0)) ? 0.0 : PetscSqr(PetscRealPart(n[5]))/PetscSqr(PetscRealPart(n[2]));
  const PetscScalar term  = 1.0 - ratio;
  return (-PetscRealPart(n[4])/(1.0 - PetscRealPart(n[3]))) - (PetscRealPart(n[2])*PetscRealPart(n[5])/(4.0*M_PI*PetscSqr(1.0 - PetscRealPart(n[3]))))*PetscSqr(term);
}


#undef __FUNCT__
#define __FUNCT__ "CalculatePhiDer_1d"
PetscErrorCode CalculatePhiDer_1d(DFTHardSphere *dft, Vec n2, Vec PhiDer2) {
  // Functions {\partial\Phi_{HS}}{\partial n_\alpha}
  typedef PetscScalar ((*PhiDerFuncType)(const PetscScalar []));
  PhiDerFuncType PhiDerFunc[6] = {PhiDer_0, PhiDer_1, PhiDer_2, PhiDer_3, PhiDer_V1, PhiDer_V2};   
  Weights *n2Array, *PhiDer2Array;
  DALocalInfo    weights2Info;
  PetscErrorCode ierr;

  PetscFunctionBegin;
  ierr = DAGetLocalInfo(dft->weightsDA2, &weights2Info);                         CHKERRQ(ierr);
  ierr = DAVecGetArray(dft->weightsDA2, n2,      &n2Array);     CHKERRQ(ierr);
  ierr = DAVecGetArray(dft->weightsDA2, PhiDer2, &PhiDer2Array);CHKERRQ(ierr);
  for(PetscInt i = weights2Info.xs; i < weights2Info.xs+weights2Info.xm; ++i) {
    for(PetscInt alpha = 0; alpha < dft->numAlpha; ++alpha) {
      PhiDer2Array[i].w[alpha] = (*(PhiDerFunc[alpha]))(n2Array[i].w);
      if (PetscIsInfOrNanScalar(PhiDer2Array[i].w[alpha])) {
        for(PetscInt beta = 0; beta < dft->numAlpha; ++beta) {
          ierr = PetscPrintf(PETSC_COMM_SELF, "alpha: %d n2Array[%d].w[%d]: %g\n", alpha, i, beta, PetscRealPart(n2Array[i].w[beta]));
        }
        SETERRQ4(PETSC_ERR_FP, "Invalid floating point value: PhiDerFunc[%d](n2[%d]) = %g+i*%g", 
                 alpha, i, PetscRealPart(PhiDer2Array[i].w[alpha]), PetscImaginaryPart(PhiDer2Array[i].w[alpha]));
      }
    }
  }
  ierr = DAVecRestoreArray(dft->weightsDA2, PhiDer2, &PhiDer2Array);CHKERRQ(ierr);
  ierr = DAVecRestoreArray(dft->weightsDA2, n2,      &n2Array);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}// CalculatePhiDer_1d()

#undef __FUNCT__
#define __FUNCT__ "DFTHardSphereCalculateMu_1d"
PetscErrorCode DFTHardSphereCalculateMu_1d(DFTHardSphere *dft, Vec rho, Vec mu){
  /*PetscErrorCode CalculateMuHS(DA da, const int numSpecies, const int numAlpha, DA singleDA, Vec PhiDer[], GeometryOptions* geomOptions, HardSphereOptions *options, Vec muHS)*/
  PetscErrorCode ierr; 
  PetscFunctionBegin;

  /* Step 1. Reflect rho.  Assume serial vector layout */
  Rho *rhoArray, *rho2Array;
  ierr = DAVecGetArray(dft->rhoDA, rho, &rhoArray); CHKERRQ(ierr);
  ierr = DAVecGetArray(dft->rhoDA2, dft->rho2, &rho2Array); CHKERRQ(ierr);
  for(PetscInt i = 0; i < dft->dims[0]; ++i) {
    for(PetscInt sp = 0; sp < dft->numSpecies; ++sp) {
      rho2Array[i].v[sp] = rhoArray[i].v[sp];
      rho2Array[dft->dims2[0]-1-i].v[sp] = rhoArray[i].v[sp];
    }
  }
  ierr = DAVecRestoreArray(dft->rhoDA, rho, &rhoArray); CHKERRQ(ierr);
  ierr = DAVecRestoreArray(dft->rhoDA2, dft->rho2, &rho2Array); CHKERRQ(ierr);

  /* Step 2. Compute n on the reflected weights DA. */
  ierr = Calculate_n_1d(dft, dft->rho2, dft->n2); CHKERRQ(ierr);

  /* Step 3. Compute PhiDer on the reflected weights DA */
  ierr = CalculatePhiDer_1d(dft, dft->n2, dft->PhiDer2); CHKERRQ(ierr);

  /* Step 4. Split PhiDer2 by alpha, put each piece in Fourier space, obtaining PhiDer2HatByAlpha.
             obtaining PhiDer2Hat.
  */
  ierr = VecStrideGatherAll(dft->PhiDer2, dft->PhiDer2ByAlpha, INSERT_VALUES); CHKERRQ(ierr);
  for(PetscInt alpha = 0; alpha < dft->numAlpha; ++alpha) { 
    ierr = MatMult(dft->F2, dft->PhiDer2ByAlpha[alpha], dft->PhiDer2HatByAlpha[alpha]);CHKERRQ(ierr); 
  }

  /* 
     Step 5. For each species, carry out the convolution of PhiDer2ByAlpha with WHat2SpeciesByAlpha,
             and sum across alphas.  Convolution is done as pointwise multiplication in Fourier space.
  */
  for(PetscInt sp = 0; sp < dft->numSpecies; ++sp) {
    ierr = VecSet(dft->mu2HatBySpecies[sp], 0.0);CHKERRQ(ierr); 
    /* Step 5a. Split WHat2BySpecies[sp] by alpha */
      ierr = VecStrideGatherAll(dft->WHat2BySpecies[sp], dft->WHat2SpeciesByAlpha, INSERT_VALUES); CHKERRQ(ierr);
      for(PetscInt alpha = 0; alpha < dft->numAlpha; ++alpha) {
        /* 
           Step 5b. Partial convolution corresponding to species sp and weight alpha. 
                    Note that z2 is a purely auxiliary vector.
         */
        ierr = VecPointwiseMult(dft->z2, dft->PhiDer2ByAlpha[alpha], dft->WHat2SpeciesByAlpha[alpha]);CHKERRQ(ierr); 
        if (dft->debug > 4) {
          PetscScalar *z2Array;
          ierr = DAVecGetArray(dft->singleDA2, dft->z2, &z2Array); CHKERRQ(ierr);
          for(PetscInt i = 0; i < dft->dims2[0]; ++i) {
            if (PetscIsInfOrNanScalar(z2Array[i])) {
              ierr = PetscPrintf(PETSC_COMM_SELF, "z2Array[%d]: %g, alpha: %d species: %d\n", i, PetscRealPart(z2Array[i]), alpha, sp);
              ierr = VecView(dft->PhiDer2ByAlpha[alpha], PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
              ierr = VecView(dft->WHat2SpeciesByAlpha[alpha], PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
              ierr = VecView(dft->z2, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
              SETERRQ(PETSC_ERR_FP, "Invalid floating point value");
            }
          }
          ierr = DAVecRestoreArray(dft->singleDA2, dft->z2, &z2Array);CHKERRQ(ierr); 
        } 
        /* Step 5c. Accumulate the piece of mu2HatBySpecies corresponding to alpha.*/
        ierr = VecAXPY(dft->mu2HatBySpecies[sp], 1.0, dft->z2);CHKERRQ(ierr); 
      }//alpha-loop
      ierr = MatMultTranspose(dft->F2, dft->mu2HatBySpecies[sp], dft->mu2BySpecies[sp]);CHKERRQ(ierr);
  }// sp-loop
  ierr = VecStrideScatterAll(dft->mu2BySpecies, dft->mu2, INSERT_VALUES);CHKERRQ(ierr);
  PetscReal scale = dft->kT/((PetscReal) dft->dims[0]);
  ierr = VecScale(dft->mu2, scale);CHKERRQ(ierr); 
  /* Step 6. Truncate mu2 --> mu.  Assume serial vector layout */
  Rho *muArray, *mu2Array;
  ierr = DAVecGetArray(dft->rhoDA, mu, &muArray); CHKERRQ(ierr);
  ierr = DAVecGetArray(dft->rhoDA2, dft->mu2, &mu2Array); CHKERRQ(ierr);
  for(PetscInt i = 0; i < dft->dims[0]; ++i) {
    for(PetscInt sp = 0; sp < dft->numSpecies; ++sp) {
      muArray[i].v[sp] = mu2Array[i].v[sp];
    }
  }
  ierr = DAVecRestoreArray(dft->rhoDA, mu, &muArray); CHKERRQ(ierr);
  ierr = DAVecRestoreArray(dft->rhoDA2, dft->mu2, &mu2Array); CHKERRQ(ierr);
  PetscFunctionReturn(0);
}// DFTHardSphereCalculateMu_1d()


#undef __FUNCT__
#define __FUNCT__ "DFTHardSphereCalculateMu"
PetscErrorCode DFTHardSphereCalculateMu(DFTHardSphere *dft, Vec rho, Vec muHS){
  PetscErrorCode ierr; 
  PetscFunctionBegin;
  if(dft->dim == 1) {
    ierr = DFTHardSphereCalculateMu_1d(dft,rho,muHS); CHKERRQ(ierr);
  }
  else {
    SETERRQ1(PETSC_ERR_SUP, "Unsupported physical domain dimension: %d", dft->dim);
  }
  PetscFunctionReturn(0);
}// DFTHardSphereCalculateMu()

